Molecular Discreteness in Reaction-Diffusion Systems Yields 
Steady States Not Seen in the Continuum Limit 
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We investigate the effects of spatial discreteness of molecules in reaction-diffusion systems. It is 
found that discreteness within the so called Kuramoto length can lead to a localization of molecules, 
resulting in novel steady states that do not exist in the continuous case. These novel states are ana- 
lyzed theoretically as the fixed points of accelerated localized reactions, an approach that was verified 
to be in good agreement with stochastic particle simulations. The relevance of this discreteness- 
induced state to biological intracellular processes is discussed. 

PACS numbers: 82.39.-k, 05.40.-a, 82.40.Ck, 87.16.-b 
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Many systems in nature that involve chemical reac- 
tions can be studied with the help of reaction-diffusion 
equations. For certain processes, a relatively small num- 
ber of suitably chosen continuous macroscopic variables 
yields excellent descriptive results. In biological systems, 
however, not only is the variety of chemicals enormous, 
the number of molecules of each of the chemical species 
can range from the relatively very large to the relatively 
very small. Now, if the species with small numbers of 
molecules were irrelevant, obviously, their existence could 
be ignored and one could focus on the species with large 
numbers of molecules that can effectively be described 
by a continuous variable. However, it should not really 
come as a surprise that it was found that, in general, 
species with small numbers of molecules cannot be ne- 
glected and that certain functions in cells can critically 
depend on very small fluctuations Q, @ . Indeed, in prior 
studies on reaction-diffusion systems some effects of fluc- 
tuations on pattern formation were found (see e.g., UQ). 
Stochastic differential equations are often used to study 
effects of fluctuations. 

Of course, on a microscopic level chemicals are com- 
posed of molecules, and the actual reactions occur be- 
tween these molecules. Therefore, in principle, reaction 
events must be integer and change only discretely. In 
analysis with stochastic differential equations, though, 
the fluctuations are regarded as continuous changes. 
Clearly, this approximation can only be valid if applied 
to fluctuations that involve sufficiently large numbers of 
molecules and should not be applied when relevant chem- 
ical species are very rare. 

In order to address this issue, we previously studied 
the effects of discreteness in simple autocatalytic reac- 
tion network systems and reported discreteness-induced 
transitions as well as drastic effects on concentrations 
A key feature of these systems was, however, that 
the medium was assumed to be well-stirred. 

In contrast, in a system with diffusion in space, the 



total number of molecules may vary from point to point. 
By assuming that the reaction is fast and the diffusion 
is slow, locally, the discreteness of the molecules can be- 
come important. In fact, this can even be the case if the 
total number of molecules is large but spread out over a 
large area as well. 

Therefore, a length scale should be considered such 
that it can serve as a benchmark for judging whether 
or not a continuum approximation is applicable. To con- 
sider this problem, the ratio between the reaction and dif- 
fusion rates is important and a candidate for the length 
scale is the typical distance over which a molecule dif- 
fuses during its lifetime, i.e., before it undergoes reaction 
as defined by Kuramoto . For reference, let us briefly 
review the work. 

Consider the reaction 
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X, 2X 



B. 



If the concentration of A is set to be constant, X is pro- 
duced at a constant rate k while decaying by the reaction 
2X — ► B at a rate k! . The average concentration of X at 
the steady state is (X) = \J kA/2k', where, for simplicity, 
A is the concentration of the chemical A. Thus the aver- 
age lifetime of X at t he stea dy state is estimated to be 
t = l/(2k'(X)) = l/V2kk'A. Suppose that X molecules 
diffuse with the diffusion constant D. The typical length 
over which an X molecule diffuses in its lifetime is then 
estimated to be 



I = V2Dt, 



(1) 
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which is called the Kuramoto length Q. 

The Kuramoto length I represents the relation between 
the reaction rate and the diffusion rate. When the system 
size is smaller than I, its behavior is dominated by diffu- 
sion and local fluctuations rapidly spread throughout the 
system. Contrastingly, if the system size is much larger 
than /, fluctuations are localized only in a small part of 
the system, and distant regions fluctuate independently. 

In this reasoning, it is assumed that the average dis- 
tance between molecules is much smaller than I. Thus 
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FIG. 1: (Color online) Time series of Ni and N 2 . r = 1, a = 4, N = 1000, L x = 1000. a) D = 10, b) D = 100, c) D = 1000. 
Initially, (jVi, N 2 , N3) = (250, 250, 500). For D = 10, X 3 reaches 0, which corresponds to the unstable fixed point (2c/3, c/3, 0). 
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FIG. 2: (Color online) Average concentration of X, for dif- 
ferent r and D (a = 4, iV = 1000, L x = 1000, sampled over 
5000 < t < 10000, and 10 trials. The error bars show the stan- 
dard deviation between the trials). The dotted lines corre- 
spond to 0.1 molecule per the Kuramoto length l\ — ^/D/50r 
for each r. 



FIG. 3: (Color online) The acceleration factor a, plotted 
against A2/J1. We measure the relation from simulations with 
different r, D, and a (N — 1000, L x — 1000, sampled over 
5000 < t < 10000, and 10 trials. The error bars show the 
standard deviation of C2 between the trials). This is very 
close to the theoretical estimation a = 1 + r-X= ■ 



the actual discreteness of the molecules can be ignored, 
and the concentration of the chemical X can be regarded 
as a continuous variable. However, if the average distance 
between molecules is comparable to or larger than I, local 
discreteness of molecules may not be negligible. Suppose 
a chemical A, with very low concentration, produces an- 
other chemical B. The average lifetime of B is short, 
such that the Kuramoto length of B is shorter than the 
average distance between adjacent A molecules. With 
this setting, chemical B molecules may be considered as 
localized around A molecules. This is especially so if the 
reactions involve 2nd or higher orders of B. Then the 
localization of chemical B may drastically alter the total 
rate of the reactions, and the effect of the local discrete- 
ness of the molecules may thus be rather significant. 

In order to systematically investigate the effects of the 
local discreteness of the molecules, we consider a simple 
one-dimensional reaction-diffusion system with 3 chemi- 



cals (Xt, X2, and X3) and the following 4 reactions 



X2 + X3 

2X 2 



>X 2 

x 2 - 



-Xi; 
X; 



x 3 - 
2X1 



-x 1 
hi 



Xi 



2X3 
-x 2 . 



Here, we assume that the first two reactions are much 
faster than the others, i.e., the reaction constants satisfy 
hi, k 2 3> k 3 > ki. To be specific, we take k\ = k 2 = lOOr, 
k 3 = ar, and fc 4 = r (r > 0, 1 < a <C 100). 

In the continuum limit, ct(t,x), the concentration of 
chemical Xj at time t and position x, is governed by the 
reaction-diffusion equation for the system given by 



dci 
~dt 
dc 2 
~di 

dc 3 
dt 



T00r(ci — C2)c3 — r{c\ — ac 2 ) 



(2) 



dx 2 



"^1 1 2 



D 



2 c 2 
dx 2 



100r(ci - c 2 )c 3 + D 3 



d 2 c 3 

dx 2 



(3) 
(4) 
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where Di is the diffusion constant of Xj,. The system is 
closed and thus the total concentration c is conserved. 
For simplicity, we assume Di = D for all i. 

The reaction-diffusion equation has fixed points at 
(ci,c 2 ,c 3 ) = (0,0,c),(VSc/(Va+ l),c/(Vo + 1),0) for 
all x. By performing a straightforward linear stability 
analysis, it is shown that only the former is stable. In- 
deed, by starting from an initial condition with c; > 0, 
this reaction-diffusion equation always converges to the 
fixed point (0, 0, c). 

Next, in order to obtain insights into the case when 
the continuum limit cannot be taken we carry out direct 
particle simulations. Each molecule diffuses randomly 
(showing Brownian motion) in a one-dimensional space 
with periodic boundary conditions (length L x ). When 
two molecules are within a distance d r they react with 
a certain probability and the total number of molecules 
(N) is conserved. 

First, we investigate the case with a = 4 and show time 
series of the number of molecules Ni of chemical species 
Xi in Fig. H As can be seen, N± and N2 do not converge 
to but to relatively large numbers. As can be expected 
the final concentrations depend on r and D and for X 2 it 
is depicted in Fig. Approximately, the concentration 
turns out to be proportional to y/r/D when N X ,N 2 *C N. 

To elucidate the origin of this proportionality, we take 
a closer look at the Kuramoto length, which, of course, 
depends on the molecule specie s. In the case of the X\ 
molecules it is given by l\ — y^D/50rc 3 , as the average 
lifetime of X\ is I/IOOC3. Here we consider the situation 
N X ,N 2 <C N, so that C3 w c. In the discussion below, we 
assume that h = y/D/btirc = y/DL x /50rN. 

Using this length l x , the density of the remaining X 2 
molecules is found to be about 0.1 molecule per l x , in- 
dependent of the parameters, as shown in Fig. [3 After 
relaxation, this density does not depend on the initial 
conditions, as long as Ni ^> 1 is satisfied initially. Fur- 
thermore, the density is independent of the system size 
L x , if L x 3> h, so that the number of remaining molecules 
N 2 is simply proportional to L x . Consequently, in this 
analysis one obtains a finite c 2 regardless of the system 
size or initial conditions which is clearly different from 
the continuum limit where c 2 goes to 0. 

In this system, X\ molecules are produced by X 2 
molecules. If A2, the average distance between X 2 
molecules, is smaller than l±, the distributions of X\ 
around neighboring X 2 molecules overlap each other sig- 
nificantly and one can regard X\ to be uniformly dis- 
tributed. In contrast, if A2 is much larger than h, 
molecules X x will localize around the X 2 molecules (The 
size L x ^> X 2 ). Then, the reaction 2Xi — > X\ + X 2 is ac- 
celerated when compared to the case that the same total 
number of X\ molecules is uniformly distributed. 

We define the acceleration factor a(A2, h) as the ratio 
between the reaction rate with localized X\ and the re- 
action rate with uniformly distributed X%. If A2 3> h, it 
is expected that a» 1. Assuming that the distribution 
of X\ is continuous and represented by the concentration 



ci(x) 13], the acceleration factor can be expressed as 
_ (4) L-ifcjdx 



(5) 



For simplicity, we assume that the distribution of the 
localized X x molecules is Gaussian with a standard de- 
viation l\ centered around the X 2 molecules (which may 
overlap each other). Suppose that the X 2 molecules are 
randomly distributed over the system with an average 
distance A2, we then obtain Q 



a = 1 



1 Aa =1 
h 



1 



20T • l x c 2 



(6) 



On the other hand, the average lifetime of X 2 
molecules is much longer, so that the Kuramoto length 
for X 2 molecules is longer than A2 . Consequently, the re- 
action 2X 2 — * X 2 +Xi is not accelerated by localization. 

Provided that N x , N 2 < ^3 , Ni w due to the fast 
reactions X 2 + X 3 -> X 2 + X x and X 3 + X x -> 2A 3 . As 
a result, the ratio between the two reaction rates is given 
by 



The rate of (X x -> X 2 ) _ ak 4 Nf 
The rate of (A 2 -> X x ) ~ 



k,N 



3^ 2 



a 
a 



(7) 



Following eq. (JTJ, the two reaction rates are balanced 
if A2 takes a value such that a — a is satisfied. Corre- 
sponding to a = a, a novel fixed point appears at 



ci = c 2 



(2(0-1)^1) X (=c s ), 



(8) 



provided c\,c 2 <C C3 and C3 = c. The stability of this 
fixed point is analyzed, by linearizing eqs. © and (JSJ 
around the fixed point. Noting that 



o 



(q-l)c. 

C2 



a- 1 



5c 2 + o{5c 2 ), (9) 



with ci — c s + Sci and c 2 = c s + Sc 2 , and rewriting eqs. 
(0) and |J2Jl with a in eq. (JjJJ, we obtain 



-2ac s - 100c (3a - l)c s + 100c 
2ac s —(3a — l)c s 

o(Sci,Sc 2 ). 



Sci 
Sc 2 

(10) 



The Jacobi matrix has two negative eigenvalues, and the 
fixed point is stable (This is natural, since if a < a, N 2 
decreases, leading to the increase of a, and vice versa). 
This fixed point (steady state) is distinct from that of 
the original reaction-diffusion equation, (0, 0, c). 

From eq. JBJ, a becomes 4 when X 2 /h — ~ 10.6. 
In our simulation with a = 4, about 0.1 X 2 molecule 
per l x remains, as shown in Fig. In other words, 
A2//1 ~ 10, in good agreement with the estimate. 

By changing a, we numerically obtain the relation be- 
tween the A2//1 and the actual acceleration factor a, 
again agreeing well with the above theoretical estimate 



a 



1 



L , as shown in Fig. |21 
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In the estimate above, we consider the case that 
Ni,N 2 <C N. On the other hand, if N is set to be smaller 
than the estimated value of N2 at the steady state, N2 
increases to satisfy the balance, and finally reaches the 
state N\ + N2 = N, A3 = 0, which corresponds to the 
unstable fixed point of the reaction-diffusion equation, 
(y/ac/(y/a + l),c/( v / a + 1),0), as shown in Fig. Q](a). 

The localization of X\ cannot be maintained without 
the spatial discreteness of X2 molecules. In reaction- 
diffusion equations, any pattern will disappear eventually 
given a sufficiently long evolution time unless it is some- 
how sustained. This is even the case when the initial 
distribution of X2 is discrete. But again, it is essential 
to recall that reaction-diffusion equations are an approx- 
imation and in that sense an idealization. In reality, a 
single molecule itself can of course not be broadened by 
diffusion and the spatial discreteness of X2 molecules is 
always maintained. By itself, a molecule is a diffusion- 
resistant pattern. 

The alteration of the steady state due to localization 
is not limited to the present type of reaction network. 
Provided that the conditions 

(i) Chemical A generates another chemical species B. 

(ii) The lifetime of B is short or the diffusion of 
B is slow so that the Kuramoto length of B is 
much smaller than the average distance between 
A molecules. 

(iii) The localization of molecule B accelerates some re- 
actions. 

are satisfied, discreteness may alter the dynamics. The 
last condition is easily satisfied if species B is involved in 
second or higher order reactions. Finally, if 

(iv) The acceleration alters the density of A molecules, 

the above acceleration mechanism may control the den- 
sity of A to produce a novel steady state. 

As for the localization effect by the discreteness of 
catalytic molecules, Shnerb et al. recently showed that 
it can amplify autocatalytic reaction-diffusion processes 



[ToL [TT| . In their model, however, the density of the cat- 
alyst is fixed as an externally given value, and the con- 
centration of the product, localized around the catalyst, 
diverges in time. In our mechanism, the density of the 
catalyst (A, or X2) changes autonomously and reaches a 
suitable value to produce the discreteness effect. Hence 
the effect of discreteness is controlled by the discreteness 
itself, leading to a novel steady state. Indeed, theoreti- 
cal estimates for the novel concentrations based on the 
self-consistent fixed point of acceleration due to the lo- 
calization agree well with numerical results. 

In so far as the conditions (i)-(iv) are met, our re- 
sult does not depend on the details of the reactions, and 
should generally be valid for reaction-diffusion systems. 
We have carried out simulations of similar reaction- 
diffusion systems, and again the discreteness effect led 
to novel pattern formation that cannot be accounted for 
by Turing type mechanisms (with or without noise). 

Experimental verification of our results should be pos- 
sible by suitably designing a reaction system, with the 
use of, say, microreactors or vesicles. Also, in biolog- 
ical cells, many chemicals work at low concentrations 
on the order of 1 nM or less. Furthermore, diffusion 
is sometimes restricted, e.g. due to surrounding macro- 
molecules, and may be slow. In such an environment, it is 
probable that the average distance between the molecules 
of a given chemical species is much larger than the Ku- 
ramoto lengths of some of the other chemical species. In- 
deed, biochemical systems contain various higher order 
reactions and positive feedback mechanisms that might 
naturally support the conditions (iii)-(iv) above. 
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